    SUBROUTINE REASHP(CSHAPE,CFLSHP,IDVOUT,IDVSHP,A1,A2,DX,X0,BETADF,PHIDF, &
                      THETADF,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP,MXN3,NAT3,  &
                      NCOMP_NEED)
      USE DDPRECISION,ONLY : WP
      IMPLICIT NONE
!                            v3
! Arguments:

      INTEGER :: IDVOUT,IDVSHP,IOSHP,MXNAT,NAT,MXN3,NAT3,NCOMP_NEED
      CHARACTER :: CDESCR*67,CFLSHP*80,CSHAPE*9
      INTEGER*2 ::    &
         ICOMP(MXNAT,3)
      INTEGER ::      &
         IXYZ(MXNAT,3)
      REAL(WP) ::        &
         A1(3),          &
         A2(3),          &
         BETADF(MXNAT),  &
         DX(3),          &
         PHIDF(MXNAT),   &
         THETADF(MXNAT), &
         X0(3)

! Local variables:

      CHARACTER :: CMSGNM*70,CWHERE*70
      CHARACTER :: DUMMY
      INTEGER :: J,K
      REAL(WP) :: R1,R2,R3,R4,R5,R6

! External routines:

      EXTERNAL ERRMSG,WRIMSG

!***********************************************************************
! Given:
!    CSHAPE = 'FROM_FILE' to read list of occupied sites and compositions
!                         from file -- use this for targets composed of
!                         isotropic materials
!             'ANIFRMFIL' to read list of occupied sites, compositions,
!                         and orientations for anisotropic material from
!                         file -- use this for targets composed of general
!                         isotropic materials
!    CFLSHP = name of (unformatted) file containing target description
!    IDVSHP = device number to use in reading this file
!    IDVOUT = device number to use for "running" output
!    IOSHP  = device number for "target.out" file
!           = -1 to suppress writing of "target.out"
!    MXNAT,MXN3 = dimensioning information
! Returns:
!    A1(3)  = target axis 1
!    A2(3)  = target axis 2
!    DX(3)  = dx/d , dy/d , dz/d   [where d=(dx*dy*dz)^{1/3}]
!    X0(3)  = location/d in TF of lattice site IXYZ=0 0 0
!    CDESCR = string describing target
!    NAT    = number of dipoles in target
!    NAT3   = 3*NAT
!    NCOMP_NEED = number of different refractive index files required
!    IXYZ(1-NAT,1-3) = lattice coordinates of dipoles
!    ICOMP(1-NAT,1-3)= composition for directions x,y,z (in Target
!                      Frame) of dipoles 1-NAT.
!    BETADF(1-NAT)=orientation angle "beta" (radians) for
!                  "Dielectric Frame" (DF) relative to Target Frame (TF)
!    PHIDF(1-NAT)=orientation angle "phi" (radians) for DF relative to T
!    THETADF(1-NAT)=orientation angle "theta" (radians) for DF relative
!                   to TF

! Note: It is expected that users may wish to modify this routine
!       to be consistent with whatever target generation software they
!       use.  However, as an example (and for test purposes) this routin
!       is supplied with code compatible with the "target.out" files whi
!       are generated by the program "calltarget"

!       Also provided below (commented out) is
!       a module of code appropriate for reading an unformatted
!       "shape.dat" file for a general target, assuming user's
!       target generation routine writes an unformatted file with
!       the same structure.


!***********************************************************************
! History:
! 91.05.02 (BTD): Added IDVOUT to argument list.
! 91.05.23 (BTD): Added A1, A2 to argument list.
! 91.11.11 (BTD): Modified for compatibility with WRITE statements in
!                 program CALLTARGET; Changed dimensioning of ICOMP from
!                 ICOMP(MXN3) to ICOMP(MXNAT,3)
! 93.09.24 (BTD): Modified to include examples of READ statements
!                 for both unformatted and formatted shape.dat files.
!                 Leave the formatted example (compatible with the
!                 "target.out" files generated by "calltarget") in effec
! 95.04.07 (BTD): Changed CDESCR*60 -> CDESCR*67 for compatibility with
!                 calling routine TARGET
! 95.05.15 (BTD): Changed  READ(DVSHP,FMT='(A60)') to ..(A67)..
!                 Added module for inhomogeneous/anisotropic targets
! 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ
! 98.02.11 (BTD): Modified to read lattice spacing DX(1-3) from
!                 "shape.dat"
! 98.03.09 (BTD): Add DX(3) to argument list
! 01.04.21 (BTD): Copied over changes made to version 5a10 on 00.11.02:
!                 Modified so now will ALWAYS expect ICOMP(x,y,z) to be
!                 included in file shape.dat
! 04.02.18 (BTD): Eliminated READ(IDVSHP,*)DX as we are not using
!                 noncubic lattices.
! 04.09.14 (BTD): Added CSHAPE,BETADF,PHIDF, and THETADF to the
!                 argument list.
!                 Add support for CSHAPE='ANIFIL' to read dipole
!                 locations, composition, and orientation angles
!                 BETADF, PHIDF, THETADF for Dielectric Frame relative
!                 to Target Frame at each occupied site.
! 07.09.11 (BTD): Changed IXYZ from INTEGER*2 to INTEGER
! 08.08.08 (BTD): Changed FRMFIL -> FROM_FILE
!                         ANIFIL -> ANIFRMFIL
! 08.08.31 (BTD): Updated for consistency with new convention for
!                 target.out files
! 08.09.12 (BTD): Removed extra read statement when called with
!                 CSHAPE=FROM_FILE
! 09.09.11 (BTD): ver7.1.0
!                 * Added NCOMP_NEED to argument list to allow checking
!                   whether enough refractive index files have been provided
! 09.09.17 (BTD): * Added string variable CWHERE
!                 * Added check to see whether size allocation is large
!                   enough for input shape file CFLSHP
!                 * Added various warning messages to tell users why
!                   failure has occurred.
! 10.01.30 (BTD): Made more robust
!                 * added diagnostics to help identify problems
!                 * read in NAT as real number in case user has
!                   written in E format
!                 * allow first "dummy" column to be character variable
!                   (DUMMY from REAL to CHARACTER)
! 10.02.07 (BTD): fixed bug
! 12.02.08 (BTD): enabled writing of target.out if IOSHP > 0
! end history

! Copyright (C) 1993,1995,1996,1998,2001,2004,2007,2008,2009,2010,2012
!               B.T. Draine and P.J. Flatau
! This code is covered by the GNU General Public License.
!***********************************************************************
!*** diagnostic
!      write(0,*)'reashp_v3 ckpt 0 mxnat=',mxnat
!***
! Enable following module of code if "shape.dat" file has been written
! in formatted form consistent with the "target.out" files generated by
! e.g., tarcel.f (inhomogeneous target, with ICOMP in file)

! Note: it is assumed that arrays THETADF,PHIDF,BETADF are
!       initialized to zero prior to calling this routine.

      CALL WRIMSG('REASHP','about to open shape file=')
      CALL WRIMSG('REASHP',CFLSHP)
      CWHERE='error opening shape file'
! diagnostic
!      write(0,*)'reashp_v3 ckpt 1 idvshp=',idvshp
!
      OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='FORMATTED',ERR=99)
! diagnostic
!      write(0,*)'reashp_v3 ckpt 2'
!
      CWHERE='error reading line 1 of shape file'
      READ(IDVSHP,FMT='(A67)',ERR=99)CDESCR                      ! read line 1
! diagnostic
!      write(0,*)'reashp_v3 ckpt 3'
!
      CWHERE='error reading NAT (line 2 of shape file)'
      READ(IDVSHP,*,ERR=99)NAT                                   ! read line 2
! diagnostic
!      write(0,*)'reashp_v3 ckpt 4, nat=',nat
!
      IF(NAT>MXNAT)THEN
! diagnostic
!         write(0,*)'reashp_v3 ckpt 4.1, mxnat=',mxnat
!
         WRITE(CMSGNM,FMT='(A,I8,A,I8,A)')                          &
              'fatal error: shape file has NAT=',NAT,               &
              ' but MXNAT=',MXNAT
         CALL WRIMSG('REASHP',CMSGNM)
         WRITE(CMSGNM,FMT='(A)')                                    &
              ' need to increase params MXNX,MXNY,MXNZ in ddscat.par'
         CALL WRIMSG('REASHP',CMSGNM)
         CALL ERRMSG('FATAL','REASHP',                               &
                     ' need to increase MNXN,MXNY,MXNZ in ddscat.par')
      ENDIF

! diagnostic
!      write(0,*)'reashp_v3 ckpt 5'
!
      CWHERE='error reading A1 (line 3 of shape file)'
      READ(IDVSHP,*,ERR=99)A1                                    ! read line 3
! diagnostic
!      write(0,*)'reashp_v3 ckpt 6'
!
      CWHERE='error reading A2 (line 4 of shape file)'
      READ(IDVSHP,*,ERR=99)A2                                    ! read line 4
! diagnostic
!      write(0,*)'reashp_v3 ckpt 7'
!
      CWHERE='error reading DX (line 5 of shape file)'
      READ(IDVSHP,*,ERR=99)DX                                    ! read line 5
! diagnostic
!      write(0,*)'reashp_v3 ckpt 8'
!
      CWHERE='error reading X0 (line 6 of shape file)'
      READ(IDVSHP,*,ERR=99)X0                                    ! read line 6
! diagnostic
!      write(0,*)'reashp_v3 ckpt 9'
!
      READ(IDVSHP,*)                                             ! read line 7
! diagnostic
!      write(0,*)'reashp_v3 ckpt 9, x0=',x0
!
      IF(CSHAPE=='FROM_FILE'.OR. &
         CSHAPE=='FRMFILPBC')THEN
! diagnostic
!         write(0,*)'reashp_v3 ckpt 10'
!
        DO J=1,NAT
            K=7+J
            WRITE(CWHERE,FMT='(A,I8,A)')               &
                 'error reading line',K,' of shape file'

            READ(IDVSHP,*,ERR=99)DUMMY,R1,R2,R3,R4,R5,R6
            IXYZ(J,1)=NINT(R1)
            IXYZ(J,2)=NINT(R2)
            IXYZ(J,3)=NINT(R3)
            ICOMP(J,1)=NINT(R4)
            ICOMP(J,2)=NINT(R5)
            ICOMP(J,3)=NINT(R6)
         ENDDO
! diagnostic
!         write(0,*)'reashp_v3 ckpt 11'
!
      ELSEIF(CSHAPE=='ANIFRMFIL'.OR. &
             CSHAPE=='ANIFILPBC')THEN
         DO J=1,NAT
            K=7+J
            WRITE(CWHERE,FMT='(A,I8,A)')               &
                 'error reading line',K,' of shape file'
            READ(IDVSHP,*,ERR=99)DUMMY,R1,R2,R3,R4,R5,R6,    &
                                 THETADF(J),PHIDF(J),BETADF(J)
            IXYZ(J,1)=NINT(R1)
            IXYZ(J,2)=NINT(R2)
            IXYZ(J,3)=NINT(R3)
            ICOMP(J,1)=NINT(R4)
            ICOMP(J,2)=NINT(R5)
            ICOMP(J,3)=NINT(R6)
         ENDDO
      ELSE
         WRITE(IDVOUT,*)'Fatal error in REASHP: ',              &
                        'do not understand target option ',CSHAPE
         STOP
      ENDIF
      CALL WRIMSG('REASHP','completed reading target file=')
      CALL WRIMSG('REASHP',CFLSHP)
      CLOSE(IDVSHP)

! check to see how many different compositions are needed
      NCOMP_NEED=1
      DO J=1,3
         DO K=1,NAT
            IF(ICOMP(K,J).GT.NCOMP_NEED)NCOMP_NEED=ICOMP(K,J)
         ENDDO
      ENDDO
! diagnostic
!      write(0,*)'reashp_v3 ckpt 20'
!      write(0,*)'    cflshp=',cflshp
!      write(0,*)'    ncomp_need=',ncomp_need
!
      WRITE(CMSGNM,FMT='(A,I3,A)')' requires', &
                       NCOMP_NEED,' different refractive index files'
      CALL WRIMSG('REASHP',CMSGNM)

! diagnostic
!      write(0,*)'reashp_v3 ckpt 21'
!
     CALL WRIMSG('REASHP',CMSGNM)

!***********************************************************************

! Enable following module of code if "shape.dat" file has been written
! in unformatted form consistent with these read statements:

!      INTEGER JX,JY,JZ
!      OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='UNFORMATTED')
!      READ(IDVSHP)CDESCR,NAT,A1,A2,DX,X0
!      READ(IDVSHP)((IXYZ(JX,JY),JX=1,NAT),JY=1,3),
!     &            ((ICOMP(J,JJ),J=1,NAT),JJ=1,3),
!     &            (THETADF(J),J=1,NAT),
!     &            (PHIDF(J),J=1,NAT),
!                  (BETADF(J),J=1,NAT)
!      CLOSE(IDVSHP)

!***********************************************************************

      NAT3=3*NAT
      IF(IOSHP>=0)THEN
         OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN')
         CWHERE='writing target.out...'
         WRITE(IOSHP,FMT=9020)NAT,A1,A2,DX,X0
         DO J=1,NAT
            WRITE(IOSHP,FMT=9030)J,IXYZ(J,1),IXYZ(J,2),IXYZ(J,3), &
                                 ICOMP(J,1),ICOMP(J,2),ICOMP(J,3)
         ENDDO
         CLOSE(UNIT=IOSHP)
      ENDIF
      RETURN
99    CALL WRIMSG('REASHP',CWHERE)
      CALL ERRMSG('FATAL','REASHP',' Error reading input shape file')
9020  FORMAT(' >REASHP: option FROM_FILE',/,                       &
         I10,' = NAT',/,                                           &
         3F10.6,' = A_1 vector',/,                                 &
         3F10.6,' = A_2 vector',/,                                 &
         3F10.6,' = lattice spacings (d_x,d_y,d_z)/d',/,           &
         3F10.5,' = lattice offset x0(1-3) = (x_TF,y_TF,z_TF)/d ', &
               'for dipole 0 0 0',/,                               &
         '     JA  IX  IY  IZ ICOMP(x,y,z)')
9030  FORMAT(I7,3I5,3I2)         
    END SUBROUTINE REASHP
